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ABSTRACT 

We have developed a 3-D Eulerian hydrodynamics code to model sources of 
gravitational radiation. The code is written in cylindrical coordinates (zcr, z, if) 
and has moving grids in the w and z-directions. We use Newtonian gravity and 
calculate the gravitational radiation in the quadrupole approximation. This 
code has been tested on a variety of problems to verify its accuracy and stability, 
and the results of these tests are reported here. 



Subject headings: hydrodynamics — methods: numerical — radiation 
mechanisms: gravitational 



1. INTRODUCTION 

Many of the most interesting astrophysical systems can be described by the equations 
of hydrodynamics coupled to gravity. Among these are sources of gravitational waves 
such as inspiralling neutron star binaries and collapsing compact stars undergoing global 
rotational instabilities (Thorne 1987, 1992; Schutz 1989). With the prospect of several 
gravitational wave detectors becoming operational within a decade (see Abramovici, et al. 
1992 and references therein), the detailed modeling of such sources has a high priority. Since 
these are time-dependent, nonlinear, and fully 3-D systems, calculating the gravitational 
radiation they produce requires numerical simulations. 

We have developed a 3-D Eulerian hydrodynamics code to model such sources. The 
code is written in cylindrical coordinates with a moving grid. It uses Newtonian gravity 
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and calculates the gravitational radiation produced in the quadrupole limit. The code 
is structured so that post-Newtonian effects, and even the full field equations of general 
relativity, can be coupled to the basic hydrodynamics routines in the future. In this paper 
we describe this code and the tests we have performed to insure its accuracy. Section |2] 
presents the basic hydrodynamical equations on which the code is built and § [3] discusses 
the numerical grids and the centering of the basic variables on those grids. The numerical 
implementation of the hydrodynamical equations is presented in detail in § |4], including a 
discussion of the finite differencing techniques used. The method used to solve Poisson's 
equation is discussed in § 0. Section |6] considers the extraction of gravitational radiation 
using various quadrupole formulae. The test-bed problems we ran on the code and the 
results we obtained are presented in § [71 Finally, § |8] summarizes our work. 



2. HYDRODYNAMICAL EQUATIONS 



We begin with the basic equations of hydrodynamics (Landau & Lifshitz 1959; Bowers 
& Wilson 1991). Conservation of mass gives the mass continuity equation 

f t+P V-v = 0, (1) 

where d/dt is the total or Lagrangian time derivative. For a moving coordinate system 

d d _ „ 

where v — v r + v g is the velocity of the fluid, v r is the velocity of the fluid relative to the 
grid, and v g is the velocity of the grid. Substituting this into ([!]) gives 

^ + V • (pv r ) + pV • v g = 0, (3) 

which is the mass continuity equation with a moving grid. The remaining equations will 
all be given in a form which shows the moving grid terms explicitly. The momentum 
conservation equation, or Euler's equation, is 

+ pvW-v g = -V-(pvv r ) - VP + /, (4) 

or, in component form, 

+ pv k V ■ v g = -V • (pu%) - V k P + f k , (5) 



- 3- 



where / is the body force density. The body force can arise from many sources, including 
self-gravity of the fluid, electric and magnetic fields generated within the fluid, and 
externally applied fields. In the current version of the code, we are only concerned with the 
body force density due to the self-gravity of the fluid. This is given by 

/ = - p v$, (6) 
where the gravitational potential $ satisfies Poisson's equation 

V 2 $ = AirGp. (7) 
We follow Bowers & Wilson (1991) and use an internal energy equation 

2^1 = -V-(peVr) - P eV-v g - PV • v, (8) 

where e is the specific internal energy of the fluid. To complete this set we need the equation 
of state (EOS) of the fluid, P = P(p,e). In this paper we use the EOS for a perfect fluid, 

P = ( 7 - l)pe. (9) 

More sophisticated equations of state, in the form of analytic expressions or tabulated 
values, can be inserted into the code in the future. 

We use a cylindrical coordinate system in our code since it gives an efficient means of 
modeling rotating systems; in addition stellar collisions and inspiralling binary stars can 
easily be set up. To formulate the hydrodynamical equations in cylindrical coordinates 
(zu,z,tp) we write v = (v,u,w), with similar expressions for v g and v r (note that w is a 
linear velocity so that w = wQ, where Q is the angular velocity). With this, and the 
definition of the gradient in cylindrical coordinates, the hydrodynamical equations take the 
following forms (Clancy 1989; Smith 1993). The mass continuity equation (^) becomes 

dp _ 1 d{wpv r ) d{pu r ) 1 d(pw r ) ( 1 d{wv g ) l du g l 1 dw g\ 
dt w dw dz w dip \zu dw dz w dtp J 

and the internal energy equation (Q) is 

d(pe) 1 d{wpev r ) d(peu r ) 1 d(pew r ) 

dt w dw dz w dtp 



w dw dz w dip J 
V w dw dz w dip I 



Euler's equation (H) in component form becomes 
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d(pv) 1 d{wpvv r ) d(pvu r ) 1 d(pvw r ) 

dt w dw dz w dip 

(\ d(wv g ) | du g | 1 _ _ 0* (12) 

\ct oct dz w dip J dw dw' 

d(pu) 1 d(wpuv r ) d(puu r ) 1 d(puw r ) 

dt w dw dz w dp 

(ld(wv g ) du g ldw g \ dP d$ 
PU ^w^w- + ^z~ + wW ~^~ P ^ (13) 



and 

d(pw) 1 d{wpwv r ) d(pwu r ) 1 d(pww r ) 

dt w dw dz w dp 

( 1 diwva) du„ 1 dw Q \ 1 dP p d$ 

V ct oct 02 w dp ) w dp w dip 

Poisson's equation (0) for the potential $ in cylindrical coordinates is 

1 d ( d$\ d 2 $ 1 <9 2 $ 

CT OCT \ OCT / OZ 2 CT 0<^ 

The EOS is independent of coordinates, and remains unchanged. Thus, (p~0|) — (P~5|) , along 
with (|9|), are the equations that will be implemented in the hydrodynamics code. 



3. NUMERICAL GRIDS AND CENTERING OF VARIABLES 



The next step is to define the numerical grids in space and time on which the 
hydrodynamical equations will be discretized. We use staggered spatial and temporal grids, 
and place the physical variables in such a way as to simplify calculations (Bowers & Wilson 
1991). 

For each of the cylindrical coordinates (ct, z, p) we define two grids: a "zone centered" 
grid, which has nodes at the centers of grid zones, and a "face centered" grid, which has 
nodes on the faces or boundaries separating grid zones. This is illustrated for the axial or 
z-grid in Figure [I]. In our notation, grid zone boundaries are denoted by integral indices 
and zone centers by half- integral indices. Thus, for an array of zone boundaries {z^}, the 
center of the kth axial zone is given by 

1 . , 

Z k+ i - ^\ Z k + Zk+l) , 
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and the grid spacings are 

A^ fc+ i = - (z k+ i - z k ) 

and 

A** = 2 - • 

The radial and angular grids are handled in a similar manner. For example, the jth radial 
zone is bounded by zuj and Wj+i, and has its center at 

1 / \ 

tJ7j + i = - [Wj + tU j+1 ) . 

The evolution of the fluid is carried out on grids that are staggered in time. For reasons 
that will become evident later, the zone centered quantities are defined at integral time 
steps, while face centered variables are defined at half-integral time steps. 

Throughout this paper, we use the subscript indices j, k, I and the superscript n to 
refer to the grids in w, z, ip, and t, respectively. Thus, a zone centered quantity S is written 

57+1^+1,,+ ! = S(w j+ i,z k+ i,<p l+ i,t n ), 
and a face centered quantity V is 

K%H = v{^z k+h ^ l+ ^). 

Any variable that has three half-integral spatial indices is defined at the center of a 3-D 
grid cell. Variables having two half-integral indices and one integral index are defined in the 
center of a zone face (the surface separating two neighboring zones). While it is occasionally 
necessary to evaluate a quantity at a zone edge (where four zones meet, specified by one 
half-integral and two integral indices), physical quantities are never evaluated at zone 
corners (all integral indices). 

In order to reduce memory requirements we have chosen to let our grid cover only the 
"top" of the coordinate system z > 0, which restricts us to treating problems for which 
z = is a plane of reflection symmetry. This is not currently a serious limitation, as there 
is a large class of interesting problems possessing such a plane of symmetry. Furthermore, 
this restriction can be removed in the future with a small number of changes. 

In our code scalar quantities, such as mass, density, and gravitational potential, are 
defined at the zone centers and at integral time steps. Vector quantities, such as velocities, 
are defined on the faces between zones and at half-integral time steps. This is the most 
natural choice for locating the variables, since the velocities are updated using the gradients 
of scalar quantities (pressure and gravitational potential), which are naturally defined on 
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the zone faces, and the velocities are used to transport scalar quantities between zones. The 
velocity components are face-centered in the coordinate along which they are directed, as 
illustrated in Figure ^. 

The physical variables thus map onto the spatial and temporal numerical grids as 
follows: 

density: p{w, z,<p,t) -> 

specific internal energy: s(w, z,<p,t) — > e n , i , i , i 
pressure: p(m,z,<p,t) -> p" + i >jt+ i^ + i 

gravitational potential: $(zn,z,(p,t) — > $ n i , i , i 

fluid velocity: v m {w, z, if, t) -> v.^x^ +1 

Th-\- — 

v g (w,z,<p,t) -> 

1 

71 -h — 



2 ,.v , 2 , 



grid velocity: v gm {w,t) -> v gj 2 ,v 



9j > ' v, . i 



2 

Vg z {Z, I) — > U gk , Ug k+ i 



Note that each component of the grid velocity maps into two 1-D arrays, one for the zone 
centered grid and one for the face centered grid in that coordinate. 

It is sometimes necessary to evaluate a zone centered quantity on a face centered grid; 
we accomplish this by using a volume average (Clancy 1989). Consider the case of averaging 
a zone centered quantity k+ i l+ i (temporal index suppressed) in the axial direction to 
obtain the value Zone centered quantities are defined to be constant throughout 

the zone, as indicated in Figure |[ The volume averaged value of q in the z-direction is 
given by 



(16) 



which reduces to a simple average if all volume elements are equal. Since the volume 
element is 

= W 3 + \ Aw 3 + \ AZk H A ^+|' ( 17 ) 

equation fll~6|) reduces to 

_ 1 g 3+ i, t -i,i + iAyi+g 3+ i t+ i i+ iA%i 

vm,k,i+h - 2 Az~ k • ( 18 > 
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Similarly, volume averaging in the angular direction gives 

i q j+ i, k+ i,i-^vi-i + g 3 -+i, fc+ i, t+ i a^ + i 

= 2 —a— a- 3 ( 19 ) 

A problem arises if this simple prescription is followed for volume averaging in the 
zu-direction. The radial averaged analog of ( |TED is 

= 2 w~Kwj ' (20) 

however, this expression will introduce inaccuracy arising from the curvature of the 
coordinate system. For example, if g J+ i = q^_i (angular and axial subscripts suppressed), 
then the face centered quantity qj should have this same value. If ( PCD is applied, however, 
we see that qj differs from the values on either side by a factor 

q. _ i {wAw)j_^ + (wAw) j+ i _ ^ ^ 1 Aw J+ i - Aw^i 



q j+ i 2 WjAzuj 2 w s 

The numerical error present in this term is particularly pronounced near the axis coordinate 
singularity. However, it can easily be avoided by replacing the radius Wj in the denominator 
of ( pUP with a curvature corrected radius Wj given by 



1 (roAro) ? . i + {wAw) j+ x 



(21) 



1 2 Aujj 
to yield 

= 2 d5~Azu~j • (22) 

Using equation ( p2[) provides a more accurate volume average in the radial direction, and 
helps to avoid singular behavior on the coordinate axis. 

It is also necessary to evaluate face centered quantities at zone centers, and this is 
accomplished using expressions analogous to (|18|) , (|19D , (|2~ID , and (p2|). These are given by 



1 g i+ i, fc+ i,;A ^ + g j+ i |H i |[+1 Ay w 



and 



= 2 ~ ^ ^]TA^ . ' (25) 
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where 

1 (zuAzu)j + (zuAzu)j + i 



= o " " A „ — • ( 26 ) 



4. FINITE DIFFERENCED HYDRO EQUATIONS 



In this section we obtain the operator split, finite difference versions of the Eulerian 
hydrodynamical equations ( [TDD " (P^D- We assume that the values of the scalar variables are 
known at time step n and that the velocities are known at time step n— 



4.1. Code Architecture 



When constructing a general purpose hydrodynamics code it is important to maintain 
flexibility and allow for future modifications by using structured programming techniques. 
Our code is divided into a large number of separate subroutines, each of which is dedicated 
to a specific task and may be replaced or modified as needed for particular astrophysical 
problems (Bowers & Wilson 1991). This modular structure allows the program to be easily 
adapted to study new situations, adding or removing different processes as desired, and 
makes it possible to change a subroutine without affecting other parts of the program. 

The central portion of our code is the "main driver" routine, which is responsible for 
determining the flow of the program during execution. It determines which subroutines 
will be included in the computation, what outputs will be required, when to write out 
data, and at what point to stop the calculation. The first task of the main driver is to call 
a setup routine that reads in various input files, sets up the coordinate grids, and either 
generates or reads in the initial data, depending on the problem being run. Once the initial 
configuration has been determined control is passed back to the main driver, which calls the 
physics driver to propagate the system forward in time. Whenever appropriate the main 
driver will call the routines that output graphics, restart dumps, and other data. 

The evolution of the physical system is carried out by the physics driver, which 
coordinates the numerical integration of the hydrodynamical equations with the calculation 
of the Newtonian gravitational field and the gravitational radiation quantities. The 
hydrodynamics computation is performed by a large number of subroutines, each devoted 
to a separate task such as the advection of a particular quantity along a given coordinate 
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direction. The gravitational potential calculation is separated into computing the boundary 
conditions and solving Poisson's equation. 

This modular structure permits the simple interchange of different methods of 
implementing certain calculations, such as differencing the advection terms, solving 
Poisson's equation, or computing the gravitational radiation. It further allows for changes 
in the physical system such as specifying a different EOS or using a 1-D or 2-D grid to 
treat problems that do not require a full 3-D calculation. One may choose to skip certain 
calculations, such as not solving Poisson's equation for the case of a pure hydrodynamics 
problem without gravity. Finally, if it becomes desirable to include additional physics, such 
as radiation transport, magneto- hydrodynamics, or the effects of nuclear processes, this can 
be done by adding new subroutines for computing these processes. 



We use the technique of operator splitting (Wilson 1979; Bowers & Wilson 1991) to 



the original set of equations and are successively applied to update the physical variables. 
Thus, the full set of equations is not solved simultaneously, and each variable is not updated 
due to the full set of physical processes in a single operation. Instead, each variable is 
partially updated under the approximation that only a single process is responsible, and 
each physical process is applied sequentially to update that variable fully. The order in 
which each term is updated can be important for the accuracy and stability of the code, 
and is generally determined by experience and experiment. The order used in our code is 
based on the accumulated experience of many people with different numerical codes, and is 
thought to be particularly good for treating shocks (Clancy 1989; Bowers & Wilson 1991). 

We first divide the solution of the hydrodynamical equations into a "Lagrangian step" , 
including Lagrangian-like processes such as acceleration and work done on a given fluid 
element, and a "transport step", which includes both transport of material through the 
grid and grid motion. The Lagrangian step is performed first, with the following sequence 
of operations: 

(a) update velocities due to pressure and gravitational forces; 

(b) update velocities and energy due to artificial viscosity, using the 
updated velocities from step (a); 

(c) update energy due to mechanical or compressional work, using the updated quantities 



4.2. Operator Splitting 



solve the hydrodynamical equations 
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from step (b). 

The transport step is then performed to update mass density, energy, and velocities due to 
the transport of fluid across zone boundaries and due to the changes in the zones arising 
from the grid motion. The transport step is broken into three parts, one for advection 
in each of the three coordinate directions. To minimize numerical artifacts, the order in 
which the transport in each direction is performed is varied with each cycle. Thus, on the 
first cycle the radial advection might be performed first, followed by the axial and angular 
advections, while on the next cycle the order might be axial-angular-radial, with the order 
changing on each successive cycle until all six permutations have been exhausted. At the 
end of the transport calculation in each direction, a momentum conservation condition is 
applied to update the velocities. For each successive substep of the transport step, the 
input variables contain the updated values from the previous substep. 



4.3. Lagrangian step 

The Lagrangian step begins by updating the terms containing the accelerations due to 
gravity and pressure gradients. For the w component of momentum, with the density held 
constant, the acceleration part of fljjD is simply 

dv dP d$ 
at aw aw 

Integrating this over the volume element AVj = WjAwjAz k+ iA(p l+ i we obtain, after some 
manipulation, 

V ■ , , l , , i — u . 



hk+±,i+\ " Awi p: 



A -in 



Awj \ J+b k +h l H *i-h k +h l +h 

where v p is used to indicate that the velocity has been partially updated to the next time 
level n + \. Here, p n , , i , , i is the radial edge centered density defined by the volume 
averaging method given by (|22]) . Similarly, (0) and (|14D become 

P n-i At n P n j+ l k+ l i + l -Pj+i.fc.i^+l 

It , i , , , i = U . , ~ ' ~ - ~ 



"i+iAH-i Az k Pj+lM+i 
Az k 



A -in 

+ — {*U^M-^+k*-ki*\) (29) 



and 
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At n ^ + i fc+ i, + i-P i+ i fc+ i,i-i 



where the z and 90 face centered densities are defined by (|TJ) and (|19D . 

At this point in the code the effects of artificial viscosity are included. The discussion 
of these terms is deferred to the next section, however, so that they may be treated in 
greater detail. 

Next we consider the Lagrangian term for mechanical or compressional work. 
Separating this term out of (|TTD yields 

%) = _ P fl^ + ^ + ijM = _, V .i (31) 
at \w ovj oz w dip J 

The finite difference version of this equation is 

A(pe) (pe) c - (pe 



At n+ 2 At n+ z 



-P n+ -*(V-v) Q , (32) 



where all quantities are evaluated at zone center, the common subscripts j + £ + § 

have been omitted for convenience, and we have used the superscript C and Q to denote 
partial updates due to compressional work and artificial viscosity, respectively. The finite 
difference velocity divergence that appears in (R21) is given by 



JTJ J "r 2 

Q _ Q Q Q 

A ^+i ™ j+ iA<p l+ i 

The pressure that appears in ([^) must be evaluated at time step n + Since pressure is 
calculated from the EOS and the density and energy are known only at time step n, we 
approximate this by 

where the partial derivatives are evaluated using the known values of p and e at time step 
n. Depending on the EOS being used, an analytic form may be available for these partial 
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derivatives, or they may be interpolated from a table. Using A(pe) = p c Ae + e n Ap, we 
convert ([32]) into an equation for the change in internal energy, 



Ae 



e c -e n 



e n Ap + [P n + \ (g) n Ap] (V ■ v) Q At n+ ^ 



(34) 



We obtain Ap and p c by using the fact that, since this is a Lagrangian process, the density 
must obey the condition 



Ap 



p c - p" 



P 



c 



-(V-v)At n+ *. 



(35) 



When these values are substituted into (|34l) along with (^) and the partial derivatives of 
the EOS, we get an explicit formula for updating the internal energy due to the effects of 
mechanical work. For the case of the ideal gas EOS (Q) we have 



e c = £ n { 1 + (V ■ v)At 



l-( 7 - 


-1) 


i + \{y -i 


)At 




1 + 1(7-1) 


\ + \{V-v)At 


(V 


• v)At 



(36) 



where the velocity divergence is evaluated using (g^ 



4.4. Artificial Viscosity 



The effects of artificial viscosity are part of the Lagrangian step, and these terms 
are actually updated in our code prior to the treatment of mechanical work as noted 
above. However, we have separated the discussion of artificial viscosity from the rest of the 
Lagrangian processes because it warrants a more thorough treatment. Artificial viscosity 
is not a true physical process, but has been developed to mimic the effects of physical 
processes in shock fronts which occur on length scales much smaller than that of the 
numerical grid. The original concept was developed by Von Neumann & Richtmeyer (1949) 
and a modern discussion can be found in Bowers & Wilson (1991). 

Artificial viscosity is designed to take the place of the microphysics of true molecular 
viscosity in shock fronts. By analogy to the equations of molecular viscosity we write the 
artificial viscosity terms in the momentum and energy equations in the form (Bowers & 
Wilson 1991): 

a -^ 1 = -MpQ" > ) < 37 > 



and 



de 
dt 



Q 



ab 



dv' 



(38) 



The quantity Q is the tensor artificial viscosity defined by 




(no summation) 



(39) 



where C q and C s h ea r are constants that specify the amount of compressional and shear 
viscosity, respectively, and are given as input to the program. We use the notation 



The definition of the diagonal terms Q aa in ( |39"1) ensures that artificial viscosity is present 
under circumstances of compression but not expansion. Currently we are not including 
shear artificial viscosity in our code, so we take C s h ea r = and limit our attention to the 
compressional terms. 

When the artificial viscosity is placed on the numerical grid the compressional terms 
are zone centered. Thus the finite difference version of the radial component of (|39|) is 



assuming that Av < 0. Note that the artificial viscosity tensor is calculated using the 
velocity that has been partially updated by the pressure and potential gradients. The other 
diagonal components, again assuming that Am < and Aw < 0, are given by 





(40) 




(41) 



and 




(42) 



With these expressions the acceleration terms 
viscosity, can easily be differenced to obtain 





(43) 
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and 



At n (pQ 



n+i 



(pQ 



w. 



(45) 



PT wAifi 

Here we have suppressed the common subscripts j + + + | wherever they appear 
and have used the superscript Q to indicate the partial update due to artificial viscosity. 



Finally, the shock heating from 



is 



At nH 



Q 



zuvj Jn -1 



AtU,:,l 
~1 



+ Q 



Q C 
^ w l+l~ w l 

wA(p l+ i 



Q z 



■ u k+i- u k 



Az 



(46) 



Again, the common subscripts J + + + | have been suppressed in each term. 
Note that although the velocities used in the shock heating have already been updated by 
(ff3|)-(H5D, the artificial viscosity tensor is the same one calculated previously and used to 
update the velocities. 



4.5. Transport Step 



The transport step includes terms involving advection and the change in zone size. It 
is divided into three parts, for transport along each of the three coordinate directions. To 
minimize numerical effects, the order in which these substeps are performed is varied on 
each cycle through all possible permutations. We use a second order monotonic advection 
scheme developed by LeBlanc (Clancy 1989; Bowers & Wilson 1991), and the consistent 
advection scheme of Norman and Winkler (1986) for angular momentum transport. 

First, consider the transport of mass density in the ^-direction. The relevant transport 
terms from ([10]) are 



dp 
dt 



d(pu r 



P- 



dUn 



dz dz 

Integrating this equation over a zone volume and rearranging terms gives 



Pk+± 



1 



At n+1 2 



U 9k+1 U 9k 



At 



i 

2 r 



Az 



fc+i 



(P*U r ) k+1 ~ (p*U r ) t 



(47) 



(48) 



where we have suppressed the indices j + |, I + |, and the superscript Az is used to indicate 
that the density has been partially updated by advection in the ^-direction. The quantity 
pi represents the density interpolated to the zone face; this interpolation is performed using 
a monotonic scheme as described in § [4.6.| below. The first term on the right hand side of 
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(f48|) can be shown to be a first order approximation to the change in density due to the 
change in the cell volume caused by the motion of the z-grid. Since the second term, which 
specifies the change in density due to the transport of matter across the zone boundaries, 
will be given to second order by the monotonic interpolation scheme specified below, the 
first order term limits the accuracy of the calculation when the grid velocity is non-zero. 
However, if the grid position at the updated time step, n+1, is known, we can replace this 
approximation by the exact expression to give 

PtU 2 = " - (P* u r)k\ . ( 49 ) 

where we have also used the updated zone spacing in the second term. Since the effects of 
the change in grid spacing are treated exactly in (49[), the accuracy is determined by the 
transport term, which is second order. This gives a higher degree of accuracy than (EM) 
when the moving grid is used (and the same accuracy when the grid is stationary), and 
therefore we use (E9|) in our code. 



The energy transport in the z-directon can be treated in similar fashion, yielding 



Az!\ i \+n+ 



(P^+i = « + f - ^TT [(P*s*u r ) k+1 - (AV) fc j , (50) 



2 



where e* is calculated using a monotonic interpolation scheme, and we have again suppressed 
the common subscripts j + |, I + |. The specific internal energy at the updated time level 
is then 

- (si) 

where p Az is defined in (^). The term (pe) c in (|50D involves the internal energy t c , which 
was updated for mechanical work in equation fl36|) in the Lagrangian step. This is given by 



{pef = P°e c , (52) 

where p c is the density updated under Lagrangian conditions for mechanical work similar 
to those used to get e c . To obtain p c we apply the Lagrangian density condition (|55]), but 
include only the effects of axial velocity; this gives 

<3 _ 



where 



Pk+h ~ i + (v.^i +l At« + i' (53) 

Q Q 

+2 Az k+l 
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and the superscript Q indicates that the velocity has been updated by artificial viscosity. 



A similar procedure is used to transport the radial and axial momentum in the 
z-direction in ( |12|) and (|i~3|) , yielding 



-1 Af+5 



k+i 



Az] 



n+l 



A 7 n+l 



(p*v*u r 



'k+1 



(p*V*U r ) f 



(54) 



where the common subscripts j, I + \ have been suppressed, and 



ur. 



.Az 



U. 



Az^ 1 



(p*u*u r 



] k+l 



(p*u*u. 



rJk- 



Pk 



(55) 



where the common subscripts j + ~, I + | have been suppressed. Several of the quantities 
in these equations are not evaluated at their usual positions on the grid, and these values 
are obtained by volume averaging as described in § |3] above. In particular, the edge 
centered velocity u r ^ k ^ l+ i in (^4|) and the zone centered velocity Utj+i^+i^+i m ( p5|) must 
be obtained from the velocity field w rj+ i fc i+ i before the monotonic scheme is applied to 
obtain v* and u* in each equation. The densities p* , , i and p* i , i . i are obtained by a 

that was calculated in 



volume average of the monotonic density p* i , . :! 

The advection of angular momentum presents special problems in finite difference 



codes, and will be discussed more fully in § [4.6.| below. At this point we simply state that 
it is advantageous to perform the advection in terms of the specific angular momentum, 
J = ww. When placed on the finite difference grid, this becomes 



J 



iV 



3zu 



(56) 



The z-transport terms from ( |T4j ) can be rewritten in terms of the specific angular momentum 

as 



d(pJ) d(pJu r ) T du g 



dt 

In finite difference form this becomes 



dz 



- pJ- 



dz 



J, 



Az 



Az 



j: 



k+\ 



At r 



k+X " k+ $ Az n+1 Az n+1 



k+i 



p*J*u, 



k+1 



(p*J* 



U, 



Pk+h 



(57) 



(58) 



where we have again suppressed the indices j + |, /. The symbol J* is used to indicate that 
we have used a "smooth" value of J in the monotonic interpolation scheme for the transport 
terms; see equation (WM below. As before, we obtain the edge centered relative velocity 
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Urj+i ki by volume averaging before applying the monotonic interpolation scheme to the 
angular quantities, and get the density p* +l kl from a volume average of the interpolated 
density in (|9|). 

During the momentum advection the density has been held constant at time level n. 
In order to ensure momentum conservation as the density is updated to time step n + 1 we 
must update the velocity by 

(59) 



v n+1 = v n 9 



,n+l 



P 

where p is evaluated at the appropriate grid location for each of the velocity components. 
In component form, this gives 



v n+l 



p 



P 



,n+l 



n+1 
1 

2 >'*>" 1 2 



U j+Lk,i+ 1 



and 



,n+l 
i 

"2 ''" 1 2 ' 



u 



p" 



,n+l 



W 



P 

n P 



(60) 



,n+l 



P 



where the face centered densities are again calculated using volume averaging. This update 
completes the transport in the axial direction. 

The advection in the radial and angular directions is accomplished in a similar fashion, 
with results as follows. 

For radial advection: 



{wAwY +1 
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At"- 



(pe 



3>1 J+2 

{wAw) n . +1 Ar +i 



{pe 



{w*p*v r ) j+1 - (ro*p*t; r ). 
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(zuAzu) n , ! 



'0 + % 



Af 



zu*p*u*v r ) j+1 - (w*p*u*v r ) 1 



(61) 

(62) 
(63) 
(64) 
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and 



jAzu 



J 



(wAw) n . + ^ At n4 

i+h l (wAw) n +\ ~ (sAro 



\n+l 



w*p*J*v 1 



J'+l 



w*p*J*v, 



J i 



(65) 



Here we have suppressed the indices fe + h, I + h throughout. The coordinate 



w* = Wj — ^v rj At 



(66) 



is the radial position at which the monotonically interpolated quantities are evaluated, as 
discussed in the next section. The velocities v r j + i, v r j ik and u rj y in flB^D, fl6"4|) and (|65|), 
respectively, are obtained by volume averaging before the monotonic interpolation scheme 
is applied, and the densities p*,i, p? k and p* , are obtained by volume averaging the density 
p* that results from monotonic interpolation in the radial direction. The specific internal 
energy e Am is obtained from (|6^ ) by using (0), with the velocity divergence term replaced 
with the radial velocity divergence 



Q C 
ZU j+1 Vj +1 - ZUjVj 



For angular advection, with the indices j + ~, k + ~ suppressed throughout: 



and 



_ n 

Pi+i ~ Ph 



At n+ 2 

n+l 

1 

2 



= if*) 



c 

'+5 A,o"+l 



At 11 ' 



l 



A^ -A^ 



(p*w r ) l+1 - {p*w r \ 



(p*e*w r ) l+1 - {p*£*w r ) l 



V4 



At r 



n+l 

1 

2 



u 



At 



■iH 



1 



l +h*Atf+l wAtf+l 



{p*v*w r ) l+1 - (p*v*w r ) l 

n 

{ P *U*W r ) l+1 - {p*U*W r ) l 
n 

Pl+h 



J 



_Acpf 



Ar H 



wAcp 



n+l 



p*J*w r 



p*j*w, 



r Jii 



Pi 



(67) 

(68) 
(69) 
(70) 

(71) 



We again evaluate w r j t i, w rlik and w rl+ i by volume averaging before applying the monotonic 
interpolation and obtain p* h p\ k and p* +l by volume averaging the density p*, which is 



monotonically interpolated in the angular direction. As before, we get e Atp from (|68| ) and 
(|53|) with the velocity divergence term now replaced with its angular analog, 



(V • v v ) 




,Q 



As noted above, the order in which the advection is carried out along each of the 
coordinate directions is varied with each time step through all six permutations. This is 
done to minimize any non-physical effects that might arise from carrying out the advection 
in the same order on each cycle. The momentum conservation calculation (|59|)-(|60|) is 
performed at the end of the transport in each direction. Caution should thus be exercised in 
interpreting each of the partially updated time levels occurring in (|48D-(|7lD. For example, 
the values of p n , e , and v® may refer to the values of those variables after advection along 
one or more directions rather than at the end of the Lagrangian step. The full time level 
update n+l indicated in fl59|) -( |6Cf) is only achieved after advection has been completed in 
all three directions. 



In transporting physical quantities from one zone to another it is necessary to obtain 
an accurate value for the quantity crossing the zone face. This is accomplished by using a 
monotonic interpolation scheme from the zone center to the zone face. The method that we 
are currently using is a second order accurate scheme developed by LeBlanc (Clancy 1989; 
Bowers & Wilson 1991), although the option currently exists in the code to use either donor 
cell advection or a monotonic scheme developed by Barton (see Centrella & Wilson 1984). 

The simplest advection scheme is the donor cell method, in which the density of the 
material that flows from one cell to another is simply taken to be the density in the cell 
that it is flowing out of. For example, if we are considering mass advection in the radial 
direction between cell j — ~ and cell j + 1, then the value of the density would be taken as 
Pj = Pj-i if v r j > and p* = p- + i if v r j < 0. The total mass transferred from cell j — \ 
to cell j + | would be Am = p*WjAzA(pv r jAt. This method is only accurate to first order 
and produces large numerical diffusion (Bowers & Wilson 1991). 

In order to achieve higher accuracy, it is preferable to interpolate zone centered 
quantities to face centers using a method that preserves monotonicity in the quantity being 
advected (van Leer 1974; Bowers & Wilson 1991). Monotonic advection prevents zones 
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from being evacuated and eliminates unstable zone to zone oscillations; see Hawley, Smarr 
and Wilson (1984b) for a study of various second order monotonic schemes. We have found 
the LeBlanc method to produce the best results in test problems. 

The LeBlanc interpolation scheme is illustrated in Figure ^ for the advection of density 
in the ^-direction; a similar treatment applies to the advection of all quantities along any 
direction. The purpose of the interpolation is to determine the value of the density at the 
zone face Zk (actually, near the zone face, as described below). We use the density in zone 
z k+ i and the two neighboring zones to get the slope of the density function, S k+ i, in that 
zone, and then use that slope to interpolate the density to the zone face. To obtain S k+ i 
we first calculate the three slopes Sl, Sr, and S 2 as shown in Figure f|. We assume that the 
densities p k _i and Pfe + | are constant in the cells k — | and k + |, respectively. The "left" 
slope Sl is then calculated by using at the zone face Zk and assuming that the density 
varies linearly from that value to the value p k+ i, as shown in Figure |j. This gives 



_ Pk+i ~ Pk-h 

L — 



A^ fc+ i/2 

A similar calculation using p k+ i at the zone face z k+1 gives the "right" slope Sr, where 



OR = 



Az k+ i/2 

The third slope, S2, is calculated by first linearly interpolating the density p fc _i from the 
zone center to the zone face Zk, and the density p fc+ | to the zone face z k +i- The slope 
between those values gives S 2 , which is 



So 



'p k+ iAz k+ z + p k+ *Az k+ i p k _iAz k+ i + p k+ iAz i 



k- 

2 " 1 2 ' '" 1 2 2 



Az k ,i \ Az k ,3 + Az k , 1 Az k ,i + Az k _i 



The slope S k+ i is then chosen to be the one from the set {Sl, Sr, £2} with the minimum 
absolute value, or zero if any two of the slopes have opposite sign. The interpolated density 
is calculated using this slope and is given by 

pl = j Pk-l + \ S *-\ ( Az ^ ~ UkAt ) Uk > . (72 ) 
p k+ i - ^S k+ i [Az k+ i + u k At) u k < 



this is the density p* k that appears in the axial mass advection equation (f49[) and the other 
axial advection equations. Note that the position at which the density is evaluated is not 
the zone edge itself, but rather a distance \u k At "upwind" from the zone edge z\.- This is 
the position of the center of the volume that will be transported across that zone boundary 
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in the time step At. To simplify this illustration, we have ignored grid motion; this may be 
included by replacing the velocity u with the relative velocity u r = u — u g in (173) and in 
Figure f|. 

A similar procedure is applied to interpolate each quantity to be advected, and the 
use of a superscript * in the transport equations of § |4.5J indicates quantities that must 
be evaluated using monotonic interpolation. Note that the coordinate w*, defined in 
equation flB6|), appears in the radial transport equations. This is the position at which the 
variable is evaluated in the monotonic interpolation scheme for transport in the Ct7-direction. 

In developing the advection equations above, we stated that it was advantageous to 
recast Euler's equations in terms of specific angular momentum, J = ww, and that the 
angular velocity components are updated using this formulation. The reason for this will 
now be explained. 

The advection of angular momentum presents special problems in finite difference 
codes. Numerical effects can cause a gradual loss of total angular momentum or a 
non-physical diffusion of angular momentum outwards from the rotation axis. These effects 
can be minimized by employing a consistent advection scheme for angular momentum 
advection (Norman & Winkler 1986; Norman, Wilson, & Barton 1980). In this method, 
we must calculate the specific angular momentum J and the angular velocity Q = w/zu, in 
addition to the linear velocity in the angular direction w, which is the variable used in the 
code. On the numerical grid J is given by (|56|) and the angular velocity is obtained from 



j 



(73) 



The key to the consistent advection method is to perform advection based on the smoothest 
of the three angular quantities, w, J and Q. Each of these quantities is interpolated to 
the zone boundary using the LeBlanc monotonic scheme described above. For the case of 
advection in the axial direction the "smooth" specific angular momentum is defined to be 
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w k zu, 
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(74) 



where we have suppressed the indices j + This smooth quantity is the one that 
appears in the angular momentum transport equations (|58D , ( |B~5D and (|7T|) . The new linear 
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velocity in the angular direction is obtained by inverting fl5Tj|) at the end of the specific 
angular momentum transport. 

The equation for obtaining the smooth specific angular momentum in the angular 
direction is completely analogous to (ff4|), with an appropriate change of indices. For 
advection in the radial direction, w* must be substituted for w in ([74]) so that the radius is 
evaluated at the same position as the interpolated quantities; see equation (|66|). 



4.7. Time Step Controls 

The finite difference equations that have been developed are used to evolve a numerical 
model of a physical system forward in time by taking discrete time steps. The time steps 
cannot be arbitrarily large, as the system of equations can become unstable (Bowers & 
Wilson 1991). We now discuss the restrictions placed on the time step to produce accurate 
and stable solutions. 

Whenever a quantity is being transported (or a signal is being propagated) through 
a grid, it is important that the time step be small enough so that the quantity or signal 
cannot cross an entire grid zone in a single time step. This is the well known Courant 
condition (Courant, Friedrichs, & Lewy 1928) given by v At/ Ax < 1 in 1-D. It leads to two 
explicit criteria that the time step in our code must satisfy. 



In an Eulerian code the fluid itself is transported through the grid, with a velocity 
v r = v — v g relative to the grid. The Courant condition for fluid transport applied in each 
of the three coordinate directions is thus 



At < (Atl) j+ i k+ i l+ i 
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At < [At] 



Az, 
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(75) 
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Note that the time step constraint varies for each zone and for each coordinate direction, so 
that ( [75]) actually specifies 3iV zones time constraints, where N zones is the total number of grid 
zones. Since all of these contraints are satisfied if we satisfy the most restrictive of them, 
we define the transport time step condition to be 

At < At T , (76) 
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where 



At 1 



mm 



1 k+y+h' 



(77) 



and the minimum is taken over the entire grid. 



Although sound waves do not physically transport fluid through the grid, they do 
propagate changes in physical variables such as density and pressure. Therefore, the 
Courant condition must be satified for the propagation of sound waves through the grid. 
The adiabatic sound speed in the fluid is 




which, for the ideal gas EOS (f|), is simply 



(7f 



(79) 



This is the velocity of the sound wave with respect to the fluid through which it travels. If 
the fluid itself is moving through the grid, then the speed of sound with respect to the grid 
varies with the direction of propagation, but has a maximum of \v r \ + c s . Thus the Courant 
condition for sound wave propagation in each of the three coordinate directions is 

Aw, 
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At < (At s z ) J+lk+hl+ , 
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(80) 
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Again, the time step is limited by the most restrictive of these conditions, so the sound 
speed Courant condition becomes 

At < At s , (81) 



where 
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(At 
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zJj+i, k+ 



ipJj+i, k+i,l-\ 



(82) 



The presence of artificial viscosity introduces an explicit diffusion into the system 
and therefore leads to an additional restriction on the time step. Considering only the 
compressional viscosity, the condition for stability in each direction is (Bowers & Wilson 
1991) 

I \ r 

(no summation). (83) 
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Satisfying this for all three coordinate directions in each grid zone again yields a set 
of SN zones constraint conditions, of which we must satisfy the most restrictive. Using 
expressions (|47]|)-(|4"2j) for the artificial viscosity tensor we define 



( ALC l)j~- - - = 



lllll] I = 



4JCV 



U, , 1 



(A^+i^+i 



(84) 
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Then the artificial viscosity constraint becomes 

At < At Q , 



where 



At Q 
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(86) 



The time step must be less than the three values defined by (|77|) , (|82|), and (pq). Note 
that the sound speed condition is always more restrictive than the transport condition 
since At s < At T , where the equality only holds in the case where pressure is independent 
of density and the sound speed c s = 0. We therefore ignore the transport constraint and 
concentrate only on the sound speed and artificial viscosity constraints (|8l|) and ([35|). 



These time step restrictions are the minimal constraint conditions based on a linearized 
stability analysis (Bowers & Wilson 1991; Press, et al. 1992). To allow for nonlinear effects 
we therefore introduce two constant parameters, C s < 1 and < 1, which are given as 
input to the program, and impose the condition 



At < min 



C s At s ,C Q At Q 



(87) 



The values chosen for the constants C s and are important. If they are too large, the 
accuracy of the simulation may be compromised; if they are too small, the number of cycles, 
and hence the CPU time, required for the calculation is increased. We usually run with 
C s w C Q sa 0.3. 

There is one final condition that should be applied. The time step is recalculated for 
each cycle, and therefore is constantly varying. If the time step changes by too much from 
one cycle to the next, irregularities can arise in the computation. To avoid this, we limit the 
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maximum change from one time step to the next by introducing another constant, C 5 > 1, 
and requiring 

At n+ 5 = min \C s At s , C Q At Q , C 5 t n ~^] . (88) 

It is usually sufficient to choose C s such that the time step cannot grow by more than 
10%-20% per cycle. 

The time step calculated in ( jS8|) is used to update zone centered quantities, defined at 
integral time levels t n . Face centered quantities are defined at half-integral time levels t n ~2 ) 
and these are updated using the time step 

At n = \ (At"-5 + Ar+s) . (89) 

When At"+2 and At n are calculated, the scalar quantities are known at time level n and 
the vector quantities are known at time level n — |, and it is these values that are used 



in equations fl80|) and fl84|) . Using conservative values for the constants in ( pq) should 
guarantee that the constraints will be satisfied for the values of the variables that are 
coincident with the time step even though they have been calculated using values that are 
up to a full time step behind. 



5. NEWTONIAN GRAVITATIONAL FIELD: POISSON'S EQUATION 



The Newtonian gravitational potential $ is obtained by solving Poisson's equation 
. Expressed in finite difference form this gives 
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Rearranging terms, we find 
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_ zu j+1 Az k+ iA(p l+ i _ WjAz k+ iA(p l+ i 

Jl A 1 J2 A 
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vj j+ iA<pi +1 ' 1 zxj j+ iAip 



Wa + 1 Am a , i Acp, , i n, + iAro. + iAw, + i 

/5 " ' h ~ ~Az k ' (92) 



and 



/7 = -(/l + /2 + / 3 + /4 + /5 + / 6 ). 

We can rewrite equation (|n]) in matrix form by relabeling the 3-D arrays Pj + i k+ i i + i 

and so that they are in 1-D form p n and $ n . For a grid with dimensions 

N mnes = N m x N z x N v , the 1-D index iV can be written 

N = L + (J - 1)N V + {K- l)N v N m . (93) 

where we have defined J = j + |, K = k + ^, and L = Z + | for convenience. With this, (j9lD 
becomes 

0N+N<pN„ + 0n-n v n^ + f 7 ®N = 4irG(pAV) N . (94) 

By inspection we see that (|94|) is in matrix form 

A A/JV <lv = 47rG(pAy) M , (95) 

where Amn is a sparse, banded matrix that is tridiagonal with 2 additional diagonal bands 
above and below the main diagonal; such a matrix is called "tridiagonal with fringes" 
(Press, et al. 1992), 

We have used the preconditioned conjugate gradient (PCG) method (Press, et al. 1992) 
to solve Poisson's equation in the form (|9~5|). We have found that this produces accurate and 
stable solutions, as will be demonstrated by the test cases reported later in this section and 
in § [f]. The PCG method was derived to solve a large, sparse system of equations in which 
the matrix Amn is symmetric and positive definite, and we have found it to perform well in 
cases where the matrix is nearly symmetric. From equations ( |9~4j ) and ( |9~2"D we see that A MN 
is exactly symmetric when a uniform grid is used, and is nearly symmetric (typically within 
a few percent) for nonuniform zoning with zoning ratios near unity. These conditions are 
always met in our simulations. 

The PCG method is an iterative technique that works as follows (Clancy 1989; Press, 
et al. 1992). The equation to be solved is written in the form 

A ■ x = b, (96) 
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where A is a matrix and the dot indicates matrix multiplication. Let 

f i = b-A-x i (97) 

be the residual, and x % be the solution on the zth iteration. Next, form 

z l = A" 1 ■ r\ (98) 

where the matrix A is called the preconditioner. Then, calculate 

P = * L C T » (99) 

p i = r + p i p i - 1 , (ioo) 

and 

a' = -. ,' ^r - (101) 

p 1 • (A • p J ) 

The improved estimate to the solution is 

x i+1 = x i + a i p\ (102) 

and the new residual is 

r i+l = f i + a(A-p i ). (103) 

To start the process, choose x to be an initial guess for the solution and use (|97|) and (j98"D 
to form f and z°, respectively. Let (3° = 0, and then proceed to iterate beginning with 
equation ( |100| ). Note that equation ( |102| ) guarantees that r t+1 in ( |103| ) is in fact the 
residual given by fl9~T|). 

We continue iterating until we have converged to the solution for $ to within a specified 
tolerance tol. We monitor 2 quantities, 



1 / ' 



(104) 



and 



E 3 » d (V 2 $ - 4ttGp) 2 

and stop iterating when either ei < tol or t2 < to/. To achieve accuracy of ~ 1% in V$, we 
need tol ^ 10~ 5 . 

There are various choices for the preconditioning matrix A. We have chosen to let A 
be the diagonal part of A, 

A = diag(A), (106) 
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to minimize the amount of storage required. This is known as diagonal scaling (Meijerink 
& Van Der Vorst 1981; Press, et al. 1992). We have implemented this method using 
only 3 additional 3-D arrays, one for each of the vectors r, z, and p; the vector A ■ p in 



equation ( 101 ) can also be stored in the array used to hold z. Two of these arrays are 
shared storage spaces with the temporary arrays needed to solve the hydrodynamical 
equations, and so solving Poisson's equation requires only one additional 3-D array. Using 
the indexing convention (|9~3"D, the coefficients (|9~2"D are calculated for each N within the loop 
to form the components of the matrices A and A as needed. Overall, this method yields a 
significant reduction in memory requirements over related methods. 

We ran comparison tests of our Poisson solver with the routine SITRSOL from the 
Cray library of scientific subroutines LIBSCI (Cray Research 1992). This routine is a 
generic solver for real, sparse linear systems of the form fl9"6|), allowing for a variety of 
solution methods and preconditioning options. We modeled uniform density dust ellipsoids 
with eccentricities up to 0.993 and compared the numerically calculated values of $ with 
the analytic expressions given by Chandrasekhar (1969); cf. § \l.3.\ below. We ran with 



32 x 32 x 32 and 64 x 32 x 64 (w, z, tp) grids, with zoning ratios in w and z varying from 1.0 
(uniformly zoned) to 1.05. The outer boundary of the grid was generally set at a distance 
of twice the largest extent of the mass distribution to produce good boundary conditions; 
see below. In terms of accuracy, all of the SITRSOL algorithms that converged produced 
results that matched the solution given by our Poisson solver. In all cases the SITRSOL 
algorithms used the full matrix Amn, without assuming that it is symmetric. Since they 
gave the same result as our routine, this is confirmation that the our Poisson solver is 
reaching the right solution even with the slight non-symmetry in Amn due to nonuniform 
zoning. When evaluated for efficiency, both in terms of memory and computational speed, 
our solver was clearly superior to the other options we investigated. Although the Cray 
routine is fully optimized, it required approximately 10 times more CPU time, regardless 
of the choice of method or the mass distribution. Furthermore, it increased the memory 
requirements of the code by a factor of 10-30, depending on the choice of method and 
preconditioner. These differences are due at least in part to the fact that our routine was 
designed for the specific problem we are solving, never requires the storage of the entire 
matrix being inverted, and takes advantage of existing temporary arrays. Since the Cray 
solver is a generic sparse matrix routine it cannot efficiently take advantage of what is 
known for the specific case at hand. 

The solution of Poisson's equation also requires specifying boundary conditions for the 
potential. These are given at the edges of the grid using a spherical multipole expansion 

n=Qm=-n ' 
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where 



Qn m = [ Y: m r n pd 3 V. 
Jv 



(108) 



Here r, 9, p are the standard spherical coordinates, the asterisk indicates the complex 
conjugate, and the integral is taken over the entire source. We carry the expansion out to 
order n = 8. 

To calculate these boundary conditions numerically, we use a set of real harmonic 
functions y nm that are related to the standard spherical harmonics by the equation 



Y nm (9, if) 



The y nm are simply the associated Legendre functions with a non-standard normalization, 



y nm {9) = (-1)' 

Thus, they obey the recursion relations 
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With this, the integral ( |108|) can be written in finite difference form as 



(109) 



where r 7 & = .Ivj 2 . i + ^ . i , = tan -1 — 2 -, and the symbol o nm indicates that we have 
used y nm instead of Y nm in (|108[) . Then (|107 ) becomes 



n=0 



m=l 



1 



ynvnXyjk) I 



;no) 
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where n max indicates the highest order multipole, currently 8 in our case, and K is used to 
indicate that we take the real part of the quantity in square brackets. Multipole moments 
q nm with n + m odd are identically zero due to the reflection symmetry through z = 0. 

Note that the boundary conditions (|110|) are determined by the source and are 
calculated in a separate subroutine prior to solving Poisson's equation. 

It is a simple matter to increase the number of terms in the expansion of the 
potential ( |110| ) at the boundaries of the grid. However, due to the discrete nature of the 
grid, we must not take the expansion out too far. This decomposition of the potential 
into multipoles is based on the orthogonality of the spherical harmonics, which must be 
preserved on the grid. Thus, the accuracy of the calculation of the boundary conditions 
depends on accurately resolving the functions used in the expansion. 

For example, consider the function cos my?, where m is equal to the number of angular 
grid zones N v . If the angular grid is equally spaced, then this function will yield a value 
of 1 at each point on the grid since the argument will always be a multiple of 2ir. For an 
axisymmetric distribution, numerically integrating this function will give nonzero values for 
multipole moments with m = N^, although all multipole moments with m^fl should be 
zero by symmetry. Therefore, we must have the order of the multipole expansion n max < N v 
at least; in general we keep n max < iV^/4. 

We expect a similar restriction on n max from the evaluation of the Legendre polynomials 
at discrete grid locations. Estimating this restriction is more difficult for two reasons. First, 
the Legendre polynomials themselves are more complicated then the simple trigonometric 
functions. Also, since the argument of the polynomials, cos 6* = z/za, is not a natural 
variable on our grid, the sampling distribution for a discrete realization of the functions is a 
complicated function of the z and w grid spacings, with some regions sampled more densely 
than others. 

We must thus take care to limit the number of terms in the multipole expansion ( fLlOj) . 
We have found that errors in the calculation of the boundary values for $ are most likely 
to occur near the z-axis. For example, for a uniform density oblate spheroid of eccentricity 
0.995 on a 64 x 32 x 32 (w, z, <p) grid with the grid boundary set at twice the equatorial 
radius, the multipole expansion (|110|) gave values for $ that matched the analytic values to 



within 1% over the entire boundary except for the zones on the z-axis. On the z-axis there 
were large errors (~ 100%) when the expansion was taken to n max = 10. These did not 
improve when additional terms in the expansion were included, but were reduced to ~ 1% 
when the order of the expansion was reduced to n max = 8. 

For a given angular resolution AL, the accuracy of the multipole expansion can be 
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improved only by moving the boundary farther from the surface of the mass distribution 
to reduce the magnitude of the higher order terms that have been discarded. We have 
found that for significantly aspherical configurations it is best to set both the w and z 
boundaries at least a factor of 2 larger than the largest dimension of the mass distribution. 
For example, when this condition was imposed for a uniform density triaxial ellipsoid with 
axis ratios 4:2:1 (in the x, y, and z-directions, respectively), the potential deviated from the 
expected value by < 1% throughout the grid on a 64 x 32 x 64 (w, z, <p) uniform grid with 
n max = 8. When the ratio of the grid boundary to the largest axis was reduced from 2 to 
1.5, maximum errors of several percent were encountered for the same parameters. 



6. EXTRACTION OF GRAVITATIONAL RADIATION 



6.1. Quadrupole Formulae 

We calculate the gravitational radiation produced in these models using the quadrupole 
approximation (Misner, Thorne and Wheeler 1973). The spacetime metric is 

9ixv = + V> ( m ) 

where /i, v — 0,1,2,3, rj^ = diag(— 1, 1, 1, 1) is the metric of flat spacetime, and 

Ifyuvl ^ 1- The gravitational wave amplitude at distance r from the source is given by the 

transverse-traceless (TT) components of the metric perturbation hij, 

tt _2G-tt 

n ab ~ r( A t ab , [ LLZ ) 

where the reduced mass quadrupole moment tensor in Cartesian coordinates is 

f ab (t) = Jp(x, t) (x a x b - \5 ab r 2 ) dV. (113) 

Here, spatial indices a,b = 1,2,3, the dot indicates a time derivative d/dt, and the integral 
is performed over the volume of the source. The TT part of the quadrupole moment is 
defined by 

•*TT •■ ^ 

I a b = PacIcdPdb — ^PabPcdlcd, (H4) 

where the projection operator P ab = 5 a b — n a rib removes the parts of the quadrupole tensor 
in the propagation direction and n a = x a /r are the components of the radial unit vector 
from the source to the field point. The gravitational wave luminosity of the source is 
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where the superscript (3) indicates 3 time derivatives, there is an implied sum on repeated 
indices, and the angle brackets are used to denote an average over several wavelengths. The 
energy emitted as gravitational radiation is 



AE = J L dt. (116) 
Finally, the angular momentum carried away by the gravitational waves is 

dJ a 2 /r(2) f (3)\ /^x 



dt "5 *^/^/' 
where e a bd is the totally antisymmetric Levi-Civita tensor. 

If we use an orthonormal spherical coordinate system defined by the unit vectors 

d ^ 1 d Id 



dr ' " r <9# ' * r sin 9 d<p ' 

with the center of mass of the source located at the origin, then the TT part of I a b has only 
four non-vanishing components. Expressed in terms of the Cartesian components ( |113j ) 
these are (Kochanek, et al. 1990) 

he = {hx cos 2 + I yy sin 2 + I xy sin 20) cos 2 6 

+I ZZ sin 2 9 — [l xz cos + I yz sin 0) sin 29, 

I<M> = I xx sin 2 + i yy cos 2 - I xy sin 20, (118) 

h<t> = I<t>e 

= \ (Iyy - I xx ) cos 9 sin 20 

+I xy cos 9 cos 20 + (I xz sin — I yz cos 0) sin 9. 

With the standard definition of the polarization basis tensors, 

e + = e e ®e e -e v ® e v , 

e x = e e <g> e v - e v <g> e e , 
the gravitational radiation waveform can then be written 

h TT = h + e + + h x e x , (119) 

where the wave amplitudes are given by 

h+ = ~(Ieo-i'<M>), (120) 

rc 4 

= (121) 

rc 4 
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All of these physical quantities associated with gravitational radiation are functions 
of at least the second time derivative of the reduced quadrupole tensor I ab . The standard 
quadrupole formula (SQF) consists of using these relations in conjunction with the definition 
( |113| ) for I ab . In an Eulerian code I ab may be calculated directly by summing over the 
grid, and the time derivatives calculated by using the finite difference approximation 
d/dt — ► A/At. However, the successive application of numerical time derivatives can 
introduce a great deal of noise into the calculated quantities, especially when the time step 
varies from cycle to cycle. To reduce this problem, Finn & Evans (1990) have developed 
two partially integrated versions of the SQF that eliminate one of the time derivatives 
by calculating I ab directly. They find that both of these methods produce much cleaner 
waveforms, with the high frequency numerical noise greatly suppressed. 



If we take the time derivative of ( 113p and use the mass continuity equation (HI) we 
obtain, after some manipulation, 

L = - j (x a x b - \r 2 5 ab ) V • (pv)dV. (122) 

This is the momentum divergence formula, which Finn & Evans (1990) refer to as QF2. 
Analytically, using ( |122| ) to calculate gravitational waveforms and luminosities is equivalent 
to using ( |113| ) directly; however, in a finite difference formulation the elimination of one 
numerical time derivative greatly increases the signal-to-noise ratio. We note that (122) 



introduces a spatial derivative, the momentum divergence, into the calculation, and that 
this derivative must be performed numerically. While this can also be expected to produce 
numerical noise, it should have little effect on the resulting value for I ab since the spatial 
derivative appears inside the volume integral and the noise will be smoothed by the 
integration process. 



The divergence in ( |122| ) can be eliminated by integrating by parts and discarding the 
surface terms at infinity. The result is 

V{ a x b ) - \5 ab (v ■ x) dV, (123) 



tab = 2 I P 

where ^ 

V( a Xb) = ^ ( V a x b + V b X a ) . 

Finn & Evans (1990) refer to this as the first moment of momentum formula, QF1. Again, 
the use of (|1 23| ) in calculating the gravitational radiation is mathematically equivalent to 
using (|113|) or ( |122j ), but in a finite difference code QF1 has several advantages. The first is 
that, like QF2, one numerical time derivative is eliminated when gravitational waveforms 
are calculated using QF1. The numerical divergence that must be calculated in QF2 has 
also been eliminated. Another advantage is that the moment arm that weights the mass 
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has been reduced from r 2 to r, which decreases the emphasis on low density material at 
large radii, while the presence of the velocity increases the emphasis on the more dynamical 
regions (Finn & Evans 1990). 



6.2. Numerical Implementation of Quadrupole Formulae 



We have implemented SQF, QF1, and QF2 in our hydrodynamics code to calculate 
the gravitational radiation. We calculate I a b (for SQF) and I a b (for QF1 and QF2) on 
each of the cycles during the run and write them to disk for additional post-processing to 
determine waveforms and luminosities. Note that on cycle n SQF produces quadrupole 
tensor components defined at time level t n+1 , while QF1 and QF2 produce time derivatives 
defined at level t n+ ^. 

The numerical integration of the SQF ( |113| ) is a simple task. We begin by defining 



x j+ i,i+i = ^. + icos^ + i, (124) 
V 3+i,i+i = tJ7 i+ isinp, + i. 

The Cartesian components of I a b are then 

f» = 2 £ P J+ i, k+ i, l+ i Uj+1,1+1 ~ l r 1+i,k+i) i , (125) 

hv = 2 E Pj+h k+i,i+i [Vj+hi+h ~ l r |+i,ft+i ) AV j+hk+i,i+i> (126) 



and 



= ^p^+y+^ (4+1 - ¥)+i^) WhmM' ^ 

hy = lyx = 2j2Pj+lk+ll+l( X V)j+l,l+l AV j+lk+l,l+l, (128) 
3'Al 



iyz = hy = 0, (129) 



where we have used the fact that the summation is over the top half of the distribution only 
(z > 0) with reflection symmetry through the x — y plane. The volume element is defined 
by dT7|) . The reduced quadrupole tensor calculated from (|125|) -( p!29| ) is mathematically 
traceless, as required. In practice, however, truncation errors can lead to a small non-zero 
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trace if all components are calculated as specified. To guarantee the tracelessness of I a b, we 
replace the calculation of I yy using ( |126| ) with 

lyy = ~ {hx + hz) ■ (130) 

The numerical implementation of QF2 presents some complications because of the 
presence of both p and v, which are defined in different positions on the grid, in the 
divergence. Recall from § |3] above that the density is zone centered, while each component 
of the velocity is face centered on a different face. In addition, zone centered quantities and 
face centered quantities are defined at different time levels. In order to obtain an accurate 
discrete representation of ( |122| ) we must perform the summation on objects that are defined 
at the same space-time points. 

At the time that we calculate the gravitational radiation, the density has been updated 
to time level n + 1 and the velocities have been updated to time level n + \. By saving the 
value of the density at the previous time step we can interpolate the density to the same 
time level as the velocity, 

•, n+ ^ _ 1 Ln i „n+l 



P j+ lk+hi+h ~ 2 [Pl+hk+y+i + p] + i,k + ii + i) ■ ( 131 ) 

Here we use the simple average because time level n + | is defined to be halfway between 
levels n and n + 1, 

With this, all of the variables needed for QF2 are available at the same time level. 

To evaluate the divergence we use the identity V • (pv) = pV -v + v- Vp. The divergence 
of v is naturally defined at the zone center, 



(V • v) i+ i k+ x l+ 1 = — I — 



J+2 1+2 
AZ ^2 

I w w j+^,k+^,i (132) 
while each component of the gradient of p is naturally face centered in that coordinate, 
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(V,p) 



and 



P j+ i, fe+ i,; + i ~ P j+ i, fe+ i,i-i 



(134) 



(135) 



The integral in ( |122| ) thus naturally breaks into four separate discrete summations over the 
numerical grid, 



r ab 



2^[(x a x b -lr 2 6 ab )p(V-v)AV 

j,k,l 

+ 2 [to - \r 2 5 ab )v{V^p) AV 

j,k,l 

+ 2 ]T [(x a x b - \r 2 5 ab )u(V z p)AV 

j,k,l 

+ 2 J] [Om - \r 2 5 ab )w{V^p)AV 

j,k,l 



(136) 



where the factor of two appears to account for integration below the x — y plane. We again 
use the symmetry of the grid to write 



^xz i zx lyz i zy 







and the tracelessness of the quadrupole tensor to obtain 

^vv i^xx ~\~ i z 



(137) 



(138) 



We therefore use 



Explicitly, these are 



to obtain only the three components I xx , I zz , and I 



j-y 



lyx- 



2E 

3,k,l 



3>k,l 

+ 2 E {(*Us - ^ r U^ ^p)*v) j+hk+hl } 



(139) 



-37- 



i« = 2 Ej(% 2 + i-h 2 + i, fc+ i)[p(v-^)A^. +ife+i/+ ,} 

j,k,l 

+ 2E{(4-|r^ i& ) [u(V z p)AV) j+l2Al+l ] 

+ 2 E {(4f| " + i) [«(V v p)AV] J+4tJH . ifI } , (140) 



and 



4, = 2E{H4, i+ i[p(V-t?)Ay] J+ i fc+ i i+ J 
+ 2 E{(^k/ + ib(V^)AF]. M) , +| } 

+ 2 EH+i^ Mv^)av]. +|iM+ J 



2 E {(^ ^MAV]^^} • (141) 



2 ' 



Here the divergence and gradient terms are obtained using ( 132|) - (|135|) . The volume 



element, AV, and the position variables, x, y, and r, are defined by ([T7|) and Q124Q , 
respectively, with appropriate indices. 

The numerical representation of QF1 is also complicated by the presence of both 
zone centered and face centered quantities in the integral. In this case, however, it is not 
possible simply to break up the summation into zone centered and face centered parts. It is 
therefore necessary to interpolate one of the quantities onto the position of the other. We 
first interpolate the density to the same time level as the velocities using ( |131| ). We then 
have the choice of evaluating the density on the zone faces, or interpolating the velocity 
components to the zone center. 

As with QF2 the symmetries of the grid and the tracelessness of I a i guarantee that 
( |137|) and ( |13S| ) still hold. The Cartesian velocity components in the x and ^/-directions are 

v x = v cos ip — w sin ip, 
v y = v simp + w cos ip. 

The zone centered representation of the components of I a b is thus 

I xx = 4y^|p zu cos (p(v cos (p — w simp) — \(vzu + uz) AV\ i i> (142) 
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I xy — 4 \ p zu cos tp(v simp + -^wcostp) — ^zuw sin 2 cp AV> 1 : 1 

j,k,l ' " J J+2' k +2' l +2 



and 



uz — ±(vw + uz) 



AV] 



(143) 



(144) 



where the zone centered velocity components are calculated by volume averaging as 
described in § [3]. A face centered discrete representation of ( |123| ) can also be obtained. In 
this case the summation again breaks up into terms centered on each face, in a manner 
similar to the separation of terms in QF2. The mass density is evaluated on each zone face 
by volume averaging, and we obtain 



4 ^2 jp \wv cos 2 ip — \zov 



j,k,l 



A V} 1 1 



A^{\pzuAV} 



j,k,l 



3+hk,l+i 



3>k,l 



(145) 



j,k,l ' 2 ' 2 

+ 4 E {l/Wcos 2 ^ - sin 2 V)}. +Kk+Kl , (146) 



and 



= -4 E {Ip^AV} + 4 £ {fpzuAV} . (147) 

We have implemented both the zone centered and the face centered methods, and found no 
noticeable difference in the results. 

The Finn & Evans (1990) formulae QF1 and QF2 produce much less noisy waveforms 
than the SQF. Nevertheless, when luminosities, which require the third derivative of 
the reduced quadrupole tensor, are calculated using these expressions we still find that 
numerical noise frequently dominates the signal. This is especially true when the time 
step is changing significantly from cycle to cycle. We have therefore implemented some 
techniques to improve the signal-to-noise ratio when the numerical derivatives are taken. 

The simplest method of taking a numerical time derivative is a simple two-point 
formula, 

, rn+l rn 
f«+2 _ 1 ab ~ 1 ab 
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When At changes a lot from cycle to cycle, this can fluctuate wildly from one time level 
to the next, causing numerical noise in the signal. These fluctuations can be reduced 
somewhat by obtaining the derivative from a higher order fit to the neighboring points. We 
have implemented a 4-point cubic fit, as follows. First define 

6t mn = t m -t r \ 

01 ab — 1 ab 1 ab' 

The time derivative can then be written 

r jran 

& = - A5fnn ~ B ( 5t ™n)\ (148) 



where 



and 



Af» ISf mn — /\f pn //itP n 

A — ' V V -B(8t mn -5t pn ), (149) 

St mn -StP n K ' y ' 

srZ n /st mn st p a l /stp" 5TZ n /st mn -8r q a 2/st'' n 



g _ st mn -stp n st mn ~sti n (150) 

~ 5tP n - 5t<i n ' 

This gives the time derivative at time level t n . The choice of the time levels m, p, 

and q in ( |1 48| ) ( P~5D| ) is essentially arbitrary, although one should ideally choose nearby 

points. To provide additional noise reduction, we perform the calculation for the four 

contiguous 4-point sets {n, m,p, q] = {n, n + 1, n + 2, n + 3}, {n, n + 1, n + 2, n — 1}, 

{n, n + 1, n — 1, n — 2}, and {n, n — 1, n — 2, n — 3}, and average the results to obtain a 

value for the derivative. This procedure is repeated for each successive time derivative. 

While using the cubic fit to obtain the derivative does provide a cleaner waveform, 
there is still a significant amount of numerical noise when several time derivatives are taken 
in succession. We have therefore found it necessary to pass the data through a filter to 
smooth it after each derivative is performed. The numerical noise in the time derivative has 
a very high frequency ~ 1 / At, whereas the signal that we are seeking is expected at much 
lower frequencies ~ 1 /to, where 

(R 3 \ 1/2 



GM 



(151) 



is the dynamical time, R is the radius, and M is the mass of the system. The relation 
t D > At is guaranteed because the Courant conditions derived in § [4.7.| require that no 



significant movement of the system can occur in a single time step. We can therefore 
smooth the data using a low pass filter with a cutoff frequency on the order of ~ 1/At, 
with little degradation of the signal. By doing this, we have been able to obtain smooth 



functions for the waveforms and luminosities; see 5 7.3. below 
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There are many methods for filtering data (see, eg., Press et al. (1992)). Because 
we are usually dealing with a large number of nonuniformly spaced data points, we use 
a "real-time" smoothing technique that consists of performing a weighted average over a 
moving window. We have made it an option to choose either a recursive (input includes 
filter output from previous cycles) or nonrecursive (input includes only unfiltered data) 
filter, although we usually use the recursive method. We commonly use a triangular 
weighting function, but other forms are available. The size of the window is variable, and 
is chosen to be at least several time steps, but smaller than any expected signal features. 
An estimate of the desired window size can easily be obtained by looking at J Q &, which has 
only small amplitude numerical noise. 



7. CODE TESTS 

We have developed a 3-D code to be used as a computational laboratory for modeling 
sources of gravitational radiation. It has been well-tested at all stages of its development 
to insure its accuracy and verify its performance (Centrella, et al. 1986). In this section we 
discuss the major test-bed problems we used and present the results of running them on 
our code. 



7.1. Riemann Shock Tube 

The Riemann shock tube is a classic problem that tests the hydrodynamics portion 
of the code, in particular the shock jump conditions, shock heating, and the resolution of 
the shock and contact discontinuities (Hawley, Smarr, & Wilson 1984a,b). Initially, the 
system consists of a fluid at rest that exists in 2 states separated by a membrane. On the 
left of the membrane is hot, dense fluid and on the right is cold, rarefied fluid. At t — the 
membrane is removed and the high density fluid flows into the low density region, creating 
a shock front moving to the right and a rarefaction wave moving to the left. 

We ran the Riemann shock tube along all 3 coordinate directions. (For the ^-direction 
we had to use a very large radius to approximate 1-D planar symmetry.) We tested several 
different advection schemes and found that the LeBlanc monotonic method discussed in 
§ [4.6.| produced the best results; for this reason we use it exclusively. Figure § shows typical 
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results from running a shock tube with 100 uniformly spaced zones in the ^-direction. The 
initial conditions for this run are (Sod 1978): 



hot, dense gas : 


P = 


1.0 g/cm 3 




e = 


2.5 erg/g 




P = 


= 1.0 erg/cm 3 


cold, rarefied gas : 


P = 


0.125 g/cm 3 




e = 


2.0 erg/g 




P = 


= 0.1 erg/cm 3 



We used the perfect fluid EOS ([]) with 7 = 1.4 and ran the shock tube in the ^-direction 
As Figure |5] shows, the code results (shown by the x symbols) accurately reproduce the 
analytic solution (shown by the solid line). Our results are essentially the same as those 
obtained by Hawley, Smarr, & Wilson (1984b) for monotonic advection. 



7.2. Spherical Polytropes 



A polytrope is a self-gravitating fluid that obeys the simple EOS 

P = kp 1 = kp 1+ ^, (152) 

where k is a constant related to the specific entropy and n is called the polytropic index. 
This can be put into the form of an ideal gas EOS (^) by defining the specific internal 
energy to be 



where the 7's that appear in (|152|) and @ are chosen to be identical. The smaller values 
of n produce "stiffer" equations of state, with n = corresponding to an incompressible 
fluid. A polytrope with index n = 1.5 (7 = 5/3) corresponds to a classical ideal gas, while 
a polytrope with n = 3 (7 = 4/3) models an ideal relativistic gas. We have used spherical 
polytropes to test the hydrodynamics coupled to Newtonian gravity. Since we have a 
cylindrical grid, this provides a non-trivial test of our code. 



7.2.1. Stability and Symmetry 
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The simplest test is to generate a stable spherical equilibrium polytrope, and verify 
that the code preserves the initial configuration. We have run this test for a 7 = 5/3 
polytrope on a 100 x 100 x 64 (zcr, z, tp) grid with zoning ratios in the range 1.01 - 1.03. 
We found that the models remained stable for > 10 6 to. We compared the gravitational 
potential calculated by our Poisson solver with the analytic values and found that the 
maximum deviation was 0.4%, just outside the surface of the mass distribution, with errors 
< 0.03% throughout more than 90% of the zones. The deviations from spherical symmetry 
were ~ 0.01%. 

We found that an initially spherical polytrope placed on our numerical grid exhibited 
small amplitude pulsations, indicating that the finite difference model was not a true 
equilibrium configuration. We therefore inserted an artificial damping term into the velocity 
equations and allowed the initial configuration to "relax" into equilibrium. The model that 
resulted when the oscillations had been damped was slightly non-spherical. The asphericity 
of the polytrope, in terms of the equatorial and polar radii, is less than a grid zone. The 
moments of inertia about the z and x-axes differed by 0.44% when the relaxation was 
performed on a 32 x 32 (w, z) 2-D axisymmetric grid. 



7.2.2. Homologous Collapse 



A polytrope with n = 3 (7 = 4/3) is marginally stable, and will collapse homologously 
if the internal pressure is slightly reduced (Goldreich & Weber 1980). We have run such 
simulations in 2-D and 3-D, with moving and stationary grids of various sizes and zoning 
ratios (Evans 1986). In all cases, the collapse remained homologous through increases in 
the central density of 10 4 -10 5 , depending on the resolution of the grid in the center of the 
polytrope. 

Figure |6] shows the results of a 2-D run with a 100 x 100 (w, z) grid and a zoning 
ratio of 1.02, and Figure [7| shows a 3-D run with a 40 x 40 x 20 (w, z, ip) grid and a zoning 
ratio of 1.04; both models use moving grids in the w and z directions. For these cases, we 
set up a spherical polytrope with mass 8.4M© and radius O.O1_R . Collapse was initiated 
by reducing pressure throughout the polytrope by 1% and imposing an inwardly directed 
homologous velocity profile v cx r with maximum velocity 10 _10 cm/sec. The moving grid is 
homologous, and is specified by choosing the grid velocity of a particular zone, ~ 1/3 of the 
way from the center to the surface, to be a fixed fraction (50% for the 2-D run and 95% for 
the 3-D run) of the fluid velocity in that zone. 
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For homologous collapse the density profile, suitably scaled, should not change 
throughout the evolution and the velocity should be purely radial, v oc r, and directed 
inward. Figure |6| shows (a) the scaled density p/p c and (6) the velocity versus the scaled 
radial distance (r/Ro)p l J z for the 2-D run. Figure (a) and (6) show the same quantities 
for the 3-D run. Here, p c is the central density and Rq is the initial polytrope radius. Each 
density graph includes plots for logp c = 9, 10, 11, 12, and 13; the 2-D plots in Figure || also 
show log p c = 14. We see that the density profiles overlay each other within the resolution of 
the graphics, except for the last trace on the graph for the 2-D collapse. From the velocity 
plot for the 2-D collapse, Figure |6] (6), we see that the core is beginning to rebound when 
p c = 10 14 g/cm 3 , so we would not expect this density profile to be the same as the other 
plots. The velocity plots show a homologous profile for the central region of the polytrope, 
including more than 95% of the mass, again except for the case p c = 10 14 g/cm 3 in the 2-D 
model. 

The sphericity of the polytrope was well preserved throughout the collapse in all runs. 
Figure |8| illustrates this for the 3-D collapse described above. Each graph displays three 
density profiles, one for density along a cut along w, one along the z-axis, and one along 
the line w — z. We see that the density profiles are indistinguishable from each other in 
both cases, demonstrating that spherical symmetry has been well-preserved. 



7.3. Homogeneous Dust Ellipsoids 

An initially uniform density, rigidly rotating, pressureless dust ellipsoid will undergo 
gravitational collapse to a pancake while maintaining a uniform density and an ellipsoidal 
figure. This process is described by a small number of ordinary differential equations 
(Miller 1974; Saenz & Shapiro 1978), which can be easily integrated numerically to obtain 
a complete description of the collapse. Since the reduced quadrupole tensor for a uniform 
ellipsoid of axes a, b, and c aligned along the x, y, and z axes, respectively, is given by the 
simple expression 

} ab = — diag(2a 2 - b s - c 2 , 2b 2 - a s - c 2 , 2c 2 - a s - b 2 ), (153) 
15 

the expected gravitational radiation from this collapse is easily obtained using the 
quadrupole formula. Finally, the gravitational potential of a uniform density ellipsoid is 
also known theoretically (Chandrasekhar 1969). The collapse of such uniform ellipsoids 
thus provides an important means of testing the hydrodynamics, gravitational radiation, 
and Poisson solver of our code. 
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Figures |9|-|10| illustrate the results of simulations of the collapse of axisymmetric dust 
ellipsoids. The plots have all been scaled to the natural scales of the problem. Thus, the 
axes a = b, and c have been scaled by the initial length of the semi- major axis a . The times 
are expressed in units of the dynamical time tp for a spherical polytrope of the same mass 
M and radius R = ao; see equation (|151[) . Finally, the gravitational radiation waveforms 
and efficiencies are scaled by a /M 2 and (M/a )~ 7 / 2 , respectively. 

Figure |9| shows the behavior of an axisymmetric oblate ellipsoid a = b uniformly 
rotating about the z-axis with Q = 0.634t £) 1 and eccentricity e = 0.1. This calculation was 
run using a 2-D 130 x 130 (w, z) grid with uniform zoning and moving grids. The evolution 
of the semi- major (equatorial) and semi-minor (polar) axes is shown in Figure |9] (a), using 
solid lines for the theoretical values and dashed lines for the code results. Figure [| (b) 
shows the evolution of the eccentricity throughout the collapse as the ellipsoid approaches a 
"pancake" with e = 1. Due to the axisymmetry, I xx = I yy = —hl zz , and all the off-diagonal 
components of I^b are zero. Thus, the gravitational radiation is all in the "+" polarization, 
and the amplitude has a simple sin 2 9 dependence. The on-axis waveform r h + {6 = ip = 0) 
is shown in Figure |] (c). Finally, Figure || (d) shows the energy emitted as gravitational 
radiation AE/Mc 2 . The gravitational radiation waveforms and efficiency dispalyed here 
were obtained using QFla. Those obtained using QFlb and QF2 are indistinquishable, 
while the waveform obtained using SQF exhibits a substantial amount of numerical noise. 
Overall, we see that the code matches the results of the theoretical calculation very well. 

In addition, we ran an axisymmetric collapse starting with eccentricity e = 0.99 and 
Q = 1.284t^ 1 using a 64 x 32 (vj,z) grid with uniform zoning and moving grids. This 
eccentricity corresponds to an axis ratio of 0.141, significantly more oblate than the axis 
ratio of 0.205 used for the bar mode instability problem of the next section. The results 
are shown in Figure [TO], which gives (a) the behavior of the axes, (b) the eccentricity, (c) 
the waveform r h + (9 — <p = 0), and (d) the energy emitted as gravitational radiation 
AE/Mc 2 . We see that the code results again agree closely with the analytic results. While 
the simulation results for the hydrodynamic quantities, such as the axes, match the analytic 
results well throughout the calculation, we can see that, towards the end of the simulation, 
when the c-axis is reduced to only a few zones in the z-direction, the gravitational wave 
calculation breaks down due to insufficient resolution. 



7.4. Bar Mode Instability 
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Rotating axisymmetric fluids are subject to global rotational instabilities that arise 
from growing modes e' im ^, where m = 2 is the so-called "bar mode". These can be 
parametrized by 

" = iwp (154) 

where T is the rotational kinetic energy and W is the gravitational potential energy; see 
Durisen & Tohline (1985) and Tassoul (1978) for reviews. The dynamical bar instability 
is driven by hydrodynamics with Newtonian gravity and develops on a timescale of about 
a rotation period. It occurs for values of (3 > (3d, where (3d ~ 0.27 for a wide range of 
rotating fluids (Durisen & Tohline 1985; Shapiro & Teukolsky 1983). The growth rate and 
eigenfrequency of the bar mode can be calculated using the linearized tensor virial equations 
(TVE) (Chandrasekhar 1969). The development of the bar instability in a differentially 
rotating polytrope with n = 3/2 for various values of (3 has been studied extensively by 
Tohline, Durisen, & McCollough (1985, hereafter TDM) using 3-D numerical simulations 
which were compared to the results of the TVE analysis. Since both their numerical and 
analytic results are available, we decided to use this case, with (3 = 0.30, as a test bed 
problem for the hydrodynamics and gravity. 

We set up the initial axisymmetric equilibrium polytrope on a cylindrical grid using the 
self-consistent field method of Smith & Centrella (1992), which is based on earlier work of 
Ostriker & Mark (1968) and Hachisu (1986). We specify a rotation law of the general form 
j = j(m(w)), where j is the specific angular momentum and m{w) is the mass interior to a 
cylinder of radius w (Ostriker & Mark 1968). For comparison with previous work we use 
the rotation law for the uniformly rotating, constant density Maclaurin spheriods, 

j(m) = |(1 - (1 - m ) 2/3 ) J/M, (155) 

where J is the total angular momentum and M is the total mass. Since polytropes do not 
have constant density, this produces differential rotation. 

The free parameters in this self-consistent field method are, in addition to the rotation 
law j(m), the polytropic index n, the central density p c , and the axis ratio R p / R, where R p 
is the polar radius and R is the equatorial radius. Upon convergence to a solution of the 
equations of hydrodynamic equilibrium, this method gives the total mass M, the polytropic 
constant k (and hence the pressure P), the angular velocity and (3. Since (3 is not 

specified initially, some experimentation with the input parameters is generally necessary 
to achieve a particular value of (3. In constructing our model with (3 ~ 0.30 we were guided 
in our choice of input parameters by the values given in Tohline, Durisen, & McCollough 
(1985). We found that using R p /R = 0.205 gives (3 = 0.301. Note that this model is highly 
flattened due to the rotation. 
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The self-consistent field code produces a rotating axisymmetric equilibrium using a 
relatively high resolution grid; this model is then interpolated onto a lower resolution grid 
for evolution by the 3-D code. To evolve the test run described here we used a 64 x 32 x 64 
(zu, z, (p) grid. Calculating T and W over this grid at the initial time gives f3 = 0.295. We 
set the grid boundaries at w — 3.85-R and z = 1.Q0R = 7.91R P . This large amount of 
initially empty space is necessary both to provide room for the star to expand as the bar 
mode grows as well as to the specify the boundary conditions for the solution of Poisson's 
equation accurately. We used a finer grid in the region initially occupied by the matter and 
a coarser grid outside. The w grid is uniform to the zone j = 30, and the z grid is uniform 
up to k = 16. The zoning is chosen such that the center of the radial zone j = 25 is at R 
and the center of the axial zone k = 9 is at R p . Outside of this uniformly zoned region, the 
zone size increases linearly with zoning ratios Awj + \/Avjj = 1.03 and Az^x/ Az^ = 1.1 
The angular grid is uniformly zoned. All grids are held fixed throughout this simulation. 
To trigger the bar instability we imposed a random perturbation with amplitude 10~ 3 on 
the density in each grid cell. 



Figure [XT] shows the development of the bar mode at 4 different times during the 
evolution of this model. The first frame, at 18. 94^, shows the model near the end of the 
linear growth regime of the bar mode. The bar structure is very apparent here, and the 
model has not yet expanded significantly beyond its initial boundaries. At 20.02to we see 
the beginning of the growth of spiral arms. By 23.26£d the model has expanded significantly 
and the spiral arms are pronounced. In the final frame, at 26.76£d, the core has returned 
to a more axisymmetric configuration and the spiral arms have expanded to about twice 
the initial radius. When evolved further, the spiral arms merge into a tenuous envelope 
surrounding the central core. These results are qualitatively similar to those obtained by 
TDM. We observe the formation of shocks during the formation and expansion of the spiral 
arms, which were not modelled in the TDM simulation. 

To study the development of the bar mode quantitataively, we analyze the density in a 
ring of fixed w and z using a complex Fourier integral 

C m {w, z) = ±- ^ p(w, <p, z^dtp, (156) 

Z7T JO 

where m = 2. The bar mode amplitude is then 

|C| = \C 2 \/C , (157) 
where C (zu, z) = p(w, z) is the mean density in the ring, and the phase angle is 

(j) = tan' 1 Re(C 2 )/Im(C 2 ). (158) 
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According to the TVE analysis, the amplitude \C\ should grow exponentially with time as 
the instability develops. 

Figure |12| shows the growth of the bar mode In \C\ versus time for the ring w = 0.362i? eg 
in zone j = 10 in the equatorial plane z — 0. We find the slope d\n \C\Jdt = 0.591^ . 
This compares favorably with the TVE result d\n.\C\/dt = 0.728 ± O.OSSt^ 1 given by 
TDM (we have converted from their units). For the eigenfrequency, which is twice the 
pattern rotation speed we find o"2 = 1.899t^ , which is very close to the TVE prediction 
o 2 = 1.892 ± 0.094^ (TDM). 

We can also compare our results with the numerical calculations of TDM. They used 
a 3-D Eulerian code written in cylindrical coordinates and required that the EOS remain 
polytropic. Their code does not solve an energy equation, and therefore has no way to 
handle self-consistently the shocks that form as the spiral arms expand. They used donor 
cell advection, which is known to be very diffusive; see § [4.6.| . And, their code assumes 
"7r-symmetry" , which means that the flow is calculated only in the angular range < (p < it 
so that only the even mode distortions are modeled. They used a 31 x 15 x 32 (w, z, (p) grid, 
with the model extending out to zone 24 in the zu-direction and zone 9 in the ^-direction, 
which is equivalent in resolution to our initial setup. They found a bar mode growth rate of 
d\n\C\/dt = 0.218 and an eigenfrequency 02 = 2.102. The substantial deviation from the 
TVE result for the growth rate is attributed to the large numerical diffusion in their code; 
see TDM for details. Overall, our code, which takes advantage of more modern techniques, 
performs significantly better on this problem. 



8. SUMMARY 



We have developed a 3-D Eulerian hydrodynamics code designed to model sources of 
gravitational radiation. This code is written in cylindrical coordinates (za,z,tp), and has 
moving grids in the w and z-directions. The hydrodynamical equations are solved using 
operator splitting (Wilson 1979; Bowers & Wilson 1991). The equations are divided into 
a "Lagrangian step", which includes Lagrangian-like processes such as accelerations and 
work done on a fluid element, and a "transport step", which includes advection of material 
through the grid and grid motion. The transport is handled using a monotonic advection 
scheme developed by LeBlanc (Clancy 1989; Bowers & Wilson 1991), with the consistent 
advection scheme for angular momentum transport (Norman & Winkler 1986; Norman, 
Wilson, & Barton 1980). 
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We use Newtonian gravity, and calculate the gravitational radiation produced using 
the quadrupole approximation. Poisson's equation is solved using an efficient and robust 
3-D matrix solver based on the preconditioned conjugate gradient method. We have 
implemented the partially integrated quadrupole formulae QF1 and QF2 developed by 
Finn & Evans (1990) and numerical filtering techniques to produce smooth gravitational 
waveforms and luminosities. 

This code has been extensively tested to verify its accuracy and stability. For example, 
the Riemann shock tube has been used to test the hydrodynamics part of the code, with 
particular attention to the shock jump conditions, shock heating, and resolution of the 
shock and contact discontinuities. The homologous collapse of spherical polytropes was 
used to test the Poisson solver and the coupling between the hydrodynamics and gravity, as 
well as the ability of the code to reproduce spherical symmetry on a cylindrical grid as the 
polytrope collapsed through many order of magnitude in density. The gravitational collapse 
of axisymmetric homogeneous dust ellipsoids expanded the tests to include the gravitational 
radiation produced. Finally, the development of the dynamical bar mode instability was 
modeled, and the results compared with theoretical calculations and previous numerical 
models. In all, our code has produced accurate and stable solution to these test-bed 
problems, and we believe it is ready to be used to model more realistic astrophysical sources 
of gravitational radiation. 

We thank R. Bowers and J. LeBlanc for many enlightening discussions and much 
valuable help. This work was supported by NSF grants PHY-8451732, PHY-8706135, 
PHY-9012383, and PHY-9208914, and by LANL. The simulations were carried out at PSC, 
NCSA, and LANL. 
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Fig. 1. — Staggered z-grid illustrating zone centered and face centered coordinate definitions. 
Zone faces are denoted by integral indices Zk and zone centers by half-integral indices z k+ i. 

Fig. 2. — Centering of variables within a grid cell. Scalar quantities S^+i are 
defined in the zone center. The velocity components v w = i^ jfc+ i jJ+ i, v z = Uj + i kl+ i, and 
v<p = lUj+i fc+i / are defined on the zone faces. (Adapted from Tohline, Durisen & McCollough 
(1985).) 

Fig. 3. — Volume averaging of zone centered variable to achieve face centered value. 
Diagonally hatched regions indicate the area for which or g J+ i fc „i is constant on 

the zone centered grid. Shaded region indicates the area on the face centered grid for which 
Qj+i,k * s t° be evaluated. The angular dimension and index have been suppressed in the 
diagram for clarity. 

Fig. 4. — Second order monotonic interpolation. The slopes Sl, Sr, and S 2 are calculated as 
described in the text. The interpolated density p* k is evaluated at a distance |-UfcAt "upwind" 
from the zone edge Zk- 

Fig. 5. — Riemann shock tube. This test was run with 100 uniformly spaced zones in the 
w-direction using the initial conditions given in the text. The solid line shows the analytic 
solution, and the x symbols show the code results. Results are given at t = 21.7 sec. All 
quantities are plotted in cgs units. 

Fig. 6. — Homologous collapse of polytrope in 2-D. A 100 x 100 (w, z) moving grid was used. 
Each graph includes plots for logp c = 9, 10, 11, 12, 13, and 14. 

Fig. 7. — Homologous collapse of polytrope in 3-D. Here, a 40 x 40 x 20 (w,z,tp) grid 
was used, with moving grids in the w and z-directions. Each graph includes plots for 
logp c = 9,10,11,12, and 13. 

Fig. 8. — Sphericity test for collapsing polytrope in 3-D. In each plot 3 density profiles are 
shown, taken along a radial (zu) cut, along the z-axis, and along the line w = z. 
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Fig. 9. — Dynamics of a collapsing, uniform density, dust ellipsoid. This figure shows an 
axisymmetric case with initial eccentricity e = 0.1 rotating around the minor (z) axis with 
angular velocity f2 = 0.634t^ 1 , run on a 2-D 130 x 130 (zu, z) grid. The solid lines show the 
theoretical values, and the dashed lines show the code results, (a) Evolution of the axes. 
(b) Eccentricity (c) Gravitational waveform, (d) Energy emitted as gravitational radiation. 



Fig. 10. — Same as Fig. || but for an initial eccentricity e = 0.99 and f2 = 1.2S4t D , run on 
a 64 x 32 (zu, z) grid. 



Fig. 11. — Development of the dynamical bar mode instability in a differentially rotating 
polytrope with (3 = 0.3. Density contours are logarithmic and are plotted at density levels 
that are 1, 2, and 3 decades down from p c . (a) t = 18.94t£i (b) t = 20.02^ (c) t = 23.2Qt D 
(d) t = 26.76tD 



Fig. 12. — The growth of the bar mode instability. The amplitude of the bar mode In \C\ at 
zu = 0.362i? e? is shown versus time. The straight line has slope d\n\C\/dt = 0.591t^ 1 . 
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